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Abstract 

This paper is concerned with the identification of the geometrical structure of the boundary 
shape for a two-dimensional boundary value problem. The output least square identifica- 
tion method is considered for estimating partially unknown boundary shapes. A numerical 
parameter estimation technique using the spline collocation method is proposed. 
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INTRODUCTION 


In the last decade, there has been much practical interest in domain identification for 
distributed parameter systems. Its application includes the optimal shape design ( See 
Pironneau [15] and the references therein ), free boundary problems arising in an oil reser- 
voir [9], non-destructive evaluations in thermal testing of materials [3], etc. In most iden- 
tification techniques, unknown domains to be identified are found by the minimization of 
output least square error functions ( OLSI - method ). Within the theoretical framework 
of optimal control, the existence of the optimal solution of the OLSI has been studied by 
Chanais [8] and Murat & Simon [14] for elliptic boundary value problems with Neumann 
and Dirichlet boundary conditions. In practical application of the domain identification, 
difficulties arise with physical domains that are often composed of non-smooth geometrical 
partsjhence the study of sensitivity analysis is of great importance ( e.g. [10] [16] [18], etc ). 
Associated with those analyses, many numerical investigations of domain identifications 
have been acomplished based on finite element methods ( See [4] [12], etc ). However, 
the computation of the derivative of the cost function with respect to the unknown do- 
main tends to be expensive and time-consuming. Such computations require intricate grid 
modification techniques for finite elements related to the decomposition of the unknown 
domain. The boundary integral equation method ( BIE ) has an advantage over the defect 
mentioned above. It does not require any decomposition of the unknown domain, only of 
its boundary. Hence the application of BIE to the domain estimation allows more sophis- 
ticated yet simpler algorithms, especially, in the case where only boundary measurement 
data are available to us. 

Our goal in domain shape identification is to develop a feasible computational method 
using the BIE method. Some previous efforts are described in [19][20] and one application 
to a wing optimization problem was illustrated in Pironneau’s book [15], Ch. 8.5. In this 
paper, we propose a computational method based on spline functions which approximate 
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T 2 (Dirichlet Condition) 


\ 



T 1 (Neumann Condition) 


Figure 1. The Spatial Domain G and the Boundary T 
the unknown boundary curve. 

To explain our approach, we restrict our attention to a 2 - D elliptic system. Let G 
( C RW ) be the open domain bounded by T. As depicted in Fig. 1, the boundary is 
decomposed into two parts, i.e., 

r = r!ur 2 . 


The system behavior on G is governed by the following Laplace equation, 

Au = 0 onG (l) 

with the mixed boundary condition 


du 

dn 

= 9i 

on r x 

(2) 

u 

= 92 

on r 2 . 

(3) 


We consider the two types of domain identification problems, 
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Case-1: Identification of Dirichlet Boundary 
Case-2: Identification of Neumann Boundary. 

The problem considered in Case-1 is to estimate the geometrical shape of the Dirichlet 
boundary T 2 from the boundary observed data on Ti, i.e., 

J/i = « ki • (4) 

In Case-2, we deal with the identification on the Neumann boundary r 2 using the boundary 
measurement on T 2 , i.e., 

» - £ • (5) 

The problem treated here involves some applications to the impedance computed tomog- 
raphy, the identification of free boundary problems, the structure design of fluid flow, 
etc. 


EXISTENCE OF SOLUTIONS 

Let 0 be a constant parametrization vector among values in a given compact set O C 
i£( n ). Throughout this paper, the parameter 0 charactrizes the unknown boundary to be 
identified. We define open sets G$ whose boundaries are given by 

Ti U r 2 (0) in Case — 1 

or 

ri(0) U r 2 in Case — 2. 

From the practical point of view, actual observed data are taken from a finite number of 
sensors allocated on the boundary. Hence we represent the corresponding output of the 
system model as follows: 
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(Case-1) 


2/i 


u lr, 

\bi h i(l?i u ) dT '’ '"'b m 


(Case-2) 


du . 

Vi ~ aS lr ’ 




where 'yr,- denotes the trace operator on T< and E* are subsets of T,-. For economy of 
notation, we rewrite (4) and (5) by 


y*iP) 


I H\U 

l 


for the Case — 1 
for the Case — 2, 


( 6 ) 


respectively. The output least square error functions for both cases are then given by 


«*(*) = £ IwW-wl (* = 1 > 2 ) 


( 7 ) 


where {y d } denotes the corresponding actual data. 


The OLSI-method is stated as follows: 

(IDP) Find the optimal parameter 6* which is the solution of 

J ( (r) = mmJ,(«) (.= 1,2). ( 8 ) 


In the sequel, we first state the existence condition of the problem for the Case-1. 


Theorem 1 ( Identification of Dirichlet boundary ) Suppose that 

(H-l ) The open sets G e depend continuously on 6 G 0 in the following sense: 
if 6 —* 0* , then the Hausdorff distance 6{Gg,G 9 ) — > 0 . 
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(H-S) There exists a bounded open set C and D such that 


C CG, CJ) 


for any 0 £ 0. 

(H-S) The sets Gg have the e -cone proprty. 

(H-4) The external input gi belongs to H?(Ti). 

(H-5) The weight functions of the measurement operator Hi satisfy 

h\ £ L~(Ei) (t = 1,2, *,m). 

Then there exists at least one solution. 

The above results can be easily obtained from Lions’ lecture notes [11] and Pironneau’s 
book [15]. 

Remark 1 The e - cone property in (H-S) means that, for Vrr £ r 2 (0) , there exists a 
direction T](x) such that, for Vz £ B(x,e) U G e , 

C(e,ri(x),z) C G e 

where C(e,rj,z) denotes the half-cone of angle e, direction rj, and vertex x intersected with 
the ball B(x, e) of center x and radius e. This property is equivalent to the corresponding 
boundary curve being Lipshitz continuous. ( See Chanais [8] ). 

For the identification in Case-2, we obtain similar results. 

Theorem 2 ( Identification of Neumann boundary ) We suppose the hypotheses (H-l) 
and (H-S) in Theorem 1 and assume that 

(H-S) ' The sets G e have C 2 -regularity property. 

(H-4) ' (fe £ t(r 2 ) 
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(H-5) • h‘eL°°(E‘) (t = 1,2, • • • ,m) 

Then there exists at least one solution. 


The proof of Theorem 2 can be obtained by extending the results of Chanais [8]. For the 
precise definition of C m - regularity property, we refer to Adams [1]. 

INTEGRAL EQUATION MODEL AND ITS NUMERICAL SCHEMES 

In this section, we replace the elliptic boundary value problem given in Section 1 by an 
equivalent integral equation on the boundary curve. A numerical method is given based 
on the spline collocation method. 

Boundary Integral Equation Model 


The Green’s representation formula yields the relation (See e. g. [6]) 


c(x 0 )«(xo) = 

Jr { u (^ -log l z_x °l) -0i lo g|z-s°| j 


dri 


+ 


/„{ 


d d\L 

1<>g ^ “ Z °^ “ dnT 1<>g \ X ~ X °l ^ dT 2 


i 


where 


c(x 0 ) = 


27T 

7T 

0 


for xo G G 

for xq E T (smooth boundary) . 
otherwise 


In the sequel, we consider the case where the boundary curve T is represented by a para- 


metric representation, 


r = 


{«(*) = (&W.6W) 


t 6 [0, 1], £ is a Jordan 
curve and |^| ± 0 


Moreover, the decomposition of the boundary T into the Dirichlet and Neumann boundary 
is assumed to be given by 


ri = {£(<) I te [o,i]> 
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r. = «W|t€[f,i]}. 

Thus, the boundary state and its flux on the boundary can be rewritten by using the arc 
length t along the curve T. Let us define 

<f>{t) = f u(£(t)) for te[0,f] 

m = ^(fW) for t e [(, 1]. 

Then the elliptic boundary value problem can be reduced to the Fredholm integral equation 
for ^ and <j> n , 

(nl - Ki)<)){s) + Li<t> n (s) = -L 2 gi{s) + K 2 g 2 {s) 

for 0 < s <t ( 9 ) 

-K 3 <}>{s) + L 3 <f> n {s) = -L 4 firi(s) + (#4 - nl)g 2 (s) 

for t < s < 1, (10) 

where 

l, = l 3 V ^‘{logifW-fWIK-) § it 

and 

l, = £{ iog\m- 

respectively. 

Numerical Scheme by Spline Collocation Method 

Many numerical methods have been proposed for the solution of integral equations, the 
most commonly used of these being Galerkin, product integration and collocation methods 
( See e. <7., [2] ). Although an asymptotic error analysis is well established for the Galerkin 
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method, its numerical implementation becomes complicated since element matrices of dis- 
cretized equations are computed through a time-consuming double integration. Owing to 
the simple numerical treatment, the collocation method for solving integral equations has 
been studied by many authors. In this paper, we use the spline collocation method for 
which the asymptotic convergence analysis is available ( See [17] ). In the sequel, we briefly 
mention this method for the models (9) and (10). 

We select an increasing sequence of mesh points 

a n,m = a n uA m 

such that 

A n = {tf = tt/N i = 0, 1, • • • , N} 

a m _ = t + i(l-t)/M i = 0,1,---,M} . 

In addition, we introduce the nodal points 

- yn with i = 

tf M = 

with j = N + 1,---,N + M. 

By S^A^) ( resp. S k (A M ) ) we denote the space of all ( k — 1) -times continuously 
differentiable splines of degree k subordinate to the partitions A N ( resp. A M ), if k > 1. 
By S 0 (A N ) ( resp. 5o(A M ) ) we denote the corresponding step functions. We approximate 
the solutions of (9) and (10) as 

<£(s) « G S k ( A n ) for 0 < s < t 

< 

. <t>n{s) « <t>n,M G S k {A M ) for t < s < 1 

Then the spline collocation method is to find 

0£a><a) € s fc (A^) x S k (A M ) 

so as to satisfy the collocation equations, 
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— -L2gi(tf’ M ) + K2g2(tf ,M ) 

for j = 1,2,---,N (11) 

= -L igi (tf’ M ) + (K, - 7rI)g 2 (tf' M ) 
for j = N + 1, • • • , N + M. 

( 12 ) 

COMPUTATIONAL METHOD FOR DOMAIN IDENTIFICATION 

In this section, using the integral equation approach as stated in the previous section, we 
consider the finite approximation for the problem ( IDP ) . 

Discretized Optimization Problem 

Let rj(t, 0) and $(t, 0) be smooth curves which specify the unknown Dirichlet and Ne umann 
boundaries, ». e., 

(Case-l) 

r 2 (0) = {rj{t,0)\t<t<i, eee 
r) € C 2 (t, 1) with 

r,(lO) = m i?(M) = £(0)}, (13) 

(Case-2) 

rip) = {£(*) = <r(M) I o < t < t, eee 

rj€zC 2 (t,l) with 

*(M) = e(i) = (i4) 
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We rewrite the system models (9) and (10) as 


where 


(Case-l) 


(Case-2) 


M&,s) 


<t> 

<t>n 


= fi(0,s) 

(* = 1 . 2 ) 


Am = 

m = 


*I-K U Ii(9) 

£•,(«) J 

+ K 2 {0)g 2 
L ~L 4 {0)g 1 + ( K 4 {0 ) - nl)g 2 J 


A 2 {0) — 

/*(«) = 


jrJ-Jfip), L x {9) 
-K 3 (0), L z 

—L 2 (0)gi + K 2 {0)g 2 
L 4 {0)gi + {K 4 - nl)g 2 J ' 


Without loss of generality, we may set observation regions as 

= {f(0 + [0,t]} 

for * = 1,2, 


( 15 ) 


(16) 


where aj and o\ are given constants. Then observation model (6) is replaced by 



W(<), 


where 



hmm 


dt 

dt 


dt, 


for A: = 1,2 * = 1,2, 


In the above equation, hj.(£(t)) denote the weighting functions of measurement opera- 
tor H k in Theorem-1 or Theorem-2 . For the state approximation, we choose the set 
of 5-splines ( See [5] ) as basis elements in 5 r (0, 1) where r is a positive integer. Let 
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{ jgAf(M) (j) }^(Af)+[r/2^ a se q Uence 0 f r _th order B-splines with the knot A N ^ M \ Then 


the solution of (15) is approximated by 


<£(s) « ^(s) = £ «f +[ (^i) |+1 5 ,>( 5 ) 


Ms) « <f>nA s ) = E 




The collocation equations for (15) yield the set of linear systems for the unknown coefficient 


vector, 


w N,M = (oc N , /3 M ) 

A?’ m (0)w n ’ m = f t N ’ M (0) ’ 

(•= 1 , 2 ) 


where A^' M and denote the corresponding element matrix and vector for the Case- 

1 and Case-2. Associated with the observation operator H,, the matrices H^' M can be 


constructed as 


I 

*' l 


Using this notation, the numerical data for the parameter 9 can be obtained by computing 

y, "•"(*) = X?- u w KM 

for i = 1,2. 

Thus we describe the computational method for the problem (IDP): 

(IDP)«" Given the data {t/d}, find the solution 0 ? G © which minimizes 

J, NM W = \ \v7 M W - »d| s (is) 

subject to (17). 


In the sequel, we discuss the computor implementation for solving the problem (IDP) M,jV . 


> 


Admissible Class of Parameters 


Solving the optimization problem (IDP) M,N , we need to specify the unknown boundary 
curve so as to satisfy the hypotheses as stated in Theorems-1 and -2. To this end, we 
approximate the unknown boundary curve by cubic B-spline functions {B i3 (t)}. In the 
sequel, we consider one simple curve as follows: 

(Case-1) 


ri(t,6) = (m {t),T)?{t,6)) 
for t < t < 1 


where 


(Case-2) 


».(<) = + 

f(M) = foW.ifM)) 

for 0 < t < t 


where 


* (<) = t-piiM°) t + e,(o) 

&M) = 

In the above equations, the Fourier coefficients {of* and \ 

through the following linear systems: 


are obtained 


Ai a M = £x 

A,/?" = 


where Ai and are given by 

r M n m 

°* 2(r^iy 

1/6, 2/3, 1/6 


1/6, 2/3, 1/6 

M rv M 

2 ( 1 - 1 ) ’ U > 2 ( 1 - 1 ) 


J 
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0.1 = [Oi,C2{t),O 2 ,---,0M^{O),O M+ i}' . 


The matrix A 2 and vector 0 2 can be described in the same way. The number of dimensions 
in 0 are thus set as 


dim{0 ) = n = 


{ 


M+l 
N + l 


for the Case — 1 
for the Case — 2 


We note that the unknown boundaries constructed by these curves satisfy the hypotheses 
(H-l) and (H-3) in Theorem-1 , and (H-l) and (H-3)' in Theorem-2. Furthermore, in order 
to assure the hypothesis (H-2) in both theorems, we impose the following constraints: 


<«"'")< ««<") (i =1,2) 


where the lower and upper bounds are given constant vectors. The precise form of the 
unknown boundary will be shown in the numerical experiments. 


Optimization Algorithm 

Under the admissible parameter class stated above, we can easily evaluate the gradient of 
the cost functional (18) with respect to 0, (i=l,2). Hence many optimization techniques are 
readily applicable to our problem. Our approach for this optimization problem is to use the 
trust region method. The trust region scheme is briefly stated as follows: Let 
be a sequence generated by this algorithm. At the current point x k , we build a model of the 
cost functional ( We usually choose a quadratic model ). Then we define a region around 
x k where we believe this model to be an adequate approximation of the functional. Using 
this model, we seek a feasible direction so as to guarantee a sufficient decrease in the model 
of cost. Once we obtain the feasible direction, the exact cost functional is evaluated at the 
new point. If its value has decreased enough, this new point is acceptable and updated as 
the next iterate. The trust region is then expanded. Otherwise the new point is rejected 
and the current trust region is reduced. The effectiveness of this algorithm is its global 
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convergence properties;namely, this algorithm makes it possible to ensure convergence to 
a critical point (optimal solution) , even from starting points ( initial guesses ) that are far 
away from the optimal solution. For detailed discussions, we refer to [13] etc. 


NUMERICAL EXPERIMENTS 


In this section, we tested our method for two examples. Although these are quite sim- 
ple, problems treated here can be easily extended to the more interesting topics, such as 
the identification of free boundary and optimal shape design problems. For the imple- 
mentation of the trust region algorithm, we used a Fortran software package created by 
Dr. R. G. Carter, ICASE ( See [7] for more details ). Test computations were carried out 
on the Gould NPl at the NASA Langley Research Center. 


Example-1 : Identification of Dirichlet Boundary 


The boundary T was decomposed into r x and r 2 by t — 0.75. The boundary inputs were 
set as 


si m) = 


i 


for t e [0.25, 0.5] 

5 


(0 for te [0., 0.5) U (0.5, 0.75] 
02 (* 7 (*)) = 0 for t G (0.75,1). 


The number of knot sequence in Example-1 was set as N = 24 and M = 8. To discretize 
the system model by the spline collocation method, we use parabolic B-splines ( i. e. r = 2 ). 
The number of sensors was taken as m = 24. The initial guesses for the parameters were 
given by 



for * = 2, ••*,8 
for * = 1,9 


Table 1 shows the estimated parameter values using the artificially generated data. Figure 2 


represents the estimated parameter function 0) and true boundary shape. 
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Example-2 : Identification of Neumann Boundary 


In this example, the boundary T was devided by t = 0.25. The boundary inputs were 
preassigned as 

0i(f(O) = O for t 6 (0., 0.25), 
f 1 for te [0.5, 0.75] 

*(€«) = 

(0 for t E [0.25, 0.5) U (0.75,1.] 

The number of knot sequence in Example-2 was set as N = 8 and M = 24. We also use 
parabolic B-spline functions for the discretized system model (17). The number of sensors 
and the initial guesses for the parameters were taken the same as in Example - 1 . Table 2 
shows the estimated parameter values and Figure 3 represents the estimated parameter 
curve $ N (t, 6) and true boundary shape. 
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Table 1. Estimated Values in Example-1 . 
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Figure 2. True Curve and Estimated Boundary in Example-1 . 
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Table 2. Estimated Values in Example- 2 . 
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Figure 3. True Curve and Estimated Boundary in Example-2 . 
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